home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / pdtoep.r < prev    next >
Text File  |  1994-12-20  |  2KB  |  54 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Symmetric positive definite Toeplitz matrix.
  4.  
  5. // Syntax:      P = pdtoep ( N , M , W , THETA ) 
  6.  
  7. // Description:
  8.  
  9. //      P is an N-by-N symmetric positive (semi-) definite (SPD)
  10. //      Toeplitz matrix, comprised of the sum of M rank 2 (or, for
  11. //      certain THETA, rank 1) SPD Toeplitz matrices. Specifically,
  12.  
  13. //            T = W(1)*T(THETA(1)) + ... + W(M)*T(THETA(M)),
  14.  
  15. //      where T(THETA(k)) has (i,j) element COS(2*PI*THETA(k)*(i-j)).
  16.  
  17. //      Defaults: M = N, W = rand(M,1), THETA = rand(M,1).
  18.  
  19. //      Reference:
  20. //       G. Cybenko and C.F. Van Loan, Computing the minimum eigenvalue of
  21. //       a symmetric positive definite Toeplitz matrix, SIAM J. Sci. Stat.
  22. //       Comput., 7 (1986), pp. 123-131.
  23.  
  24. //    This file is a translation of pdtoep.m from version 2.0 of
  25. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  26. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  27.  
  28. //-------------------------------------------------------------------//
  29.  
  30. pdtoep = function ( n, m, w, theta )
  31. {
  32.   local ( n, m, w, theta )
  33.   global (pi)
  34.  
  35.   if (!exist (m)) { m = n; }
  36.   if (!exist (w)) { w = rand (m, 1); }
  37.   if (!exist (theta)) { theta = rand (m, 1); }
  38.  
  39.   if (max (size (w)) != m || max (size (theta)) != m)
  40.   {
  41.     error ("Arguments W and THETA must be vectors of length M.");
  42.   }
  43.  
  44.   T = zeros(n,n);
  45.   E = 2*pi*( (1:n)'*ones(1,n) - ones(n,1)*(1:n) );
  46.  
  47.   for (i in 1:m)
  48.   {
  49.     T = T + w[i] * cos( theta[i]*E );
  50.   }
  51.  
  52.   return T;
  53. };
  54.